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Abstract 

We present a new splitting method for time-dependent convection-dominated diffusion 
problems. The original convection diffusion system is split into two sub-systems: a pure 
convection system and a diffusion system. At each time step, a convection problem and a 
diffusion problem are solved successively. The scheme has the following nice features: the 
convection subproblem is solved explicitly and a multistep technique is introduced to essen- 
tially enlarge the stability region so that the resulting scheme behaves like an unconditionally 
stable scheme; the diffusion subproblem is always self-adjoint and coercive so that it can be 
solved efficiently using many existing optimal preconditioned iterative solvers. The scheme 
is then extended for Navicr-Stokcs equations, where the nonlinear convection is resolved by 
a linear explicit multistep scheme at the convection step, and only a generalized Stokes prob- 
lem is needed to solve at the diffusion step with the resulting stiffness matrix being invariant 
in the time marching process. The new schemes are all free from tuning some stabiliza- 
tion parameters for the convection-dominated diffusion problems. Numerical simulations are 
presented to demonstrate the stability, convergence and performance of the single-step and 
multistep variants of the new scheme. 

Key Words. Convection-dominated diffusion problems, Navier-Stokes equations, operator 
splitting, finite elements, multistep scheme. 

AMS Classification. 65M12, 65M60, 76D05 

1 Introduction 

In this work we shall propose a new fully discrete splitting scheme for solving the convection- 
dominated diffusion problems of the following general form 



ut + V ■ (bn) - V • (eVii) + cu = F in 17 x (0, T) (1) 
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with the boundary and initial conditions 

u = Ub on 9$7x(0,r); n(0, x) = tio(x) in $7 (2) 

where Q is an open bounded polyhedral domain in M"^ {d = 1, 2, 3) with boundary F = 50, and 
T is the terminal time. Functions b, e and c in ([T]) are the convective field, diffusion and reactive 
coefficients respectively, while F, Ub and uq are the source term, the boundary and initial data 
respectively. As we are mainly interested in the construction of numerical schemes, we will not 
specify detailed regularity conditions on all these coefficients to ensure the well-posedness of the 
initial-boundary value problem ([I])-(l2|). 

The new fully discrete splitting scheme is then extended for Navier-Stokes equations 

ut + (u • V)u - iie"^ Au + Vp = F in $7x(0,r) 

V-u = in nx{0,T) 

with the boundary and initial conditions 

u = Uf, on dQx{0,T); u(0,x)=uo(x) in (4) 

where u, p, F and Re are respectively the velocity, pressure, body force and Reynolds number, 
while Ufe and Uq are the given boundary and initial data. 

The numerical solution of a time-dependent problem requires a discretization in both time 
and space, and some linearization if the concerned problem is nonlinear. A great variety of time 
marching schemes are available in the literature, such as the classical methods like the forward 
and backward Euler schemes, the Crank-Nicolson scheme, the Adams-Bashforth method etc. 
Operator splitting is also a popular technique for time discretization, such as the Yanenko 
method, the Peaceman-Rachford method, the Douglas-Rachford method and the 6 scheme; see 
[H [21 [3] and references therein. 

In solving the convection-dominated diffusion equations and the Navier-Stokes equations 
with large Reynolds numbers, it is known that standard finite element methods perform poorly 
and may exhibit nonphysical oscillations. Many spatial stabilization techniques have been pro- 
posed and studied. The streamline-upwind Petrov-Galerkin method was developed for convective 
transport problems [11[5], and its basic idea is to modify the standard Petrov-Galerkin formu- 
lation by adding a streamline upwind perturbation, which acts only in the flow direction and is 
solely defined in the interiors of elements. The Galerkin least-squares method [6] is a conceptual 
simplification of the streamline-upwind Petrov-Galerkin method, and adds a stabilization that 
involves an element-by-element weighted least-squares of the residual to the original differential 
equation. The efficiency of these two stabilization techniques may be affected by the choices of 
the stabilization parameters involved. There are still no precise general formulae to help select 
optimal parameters in numerical simulations; see, e.g., [Tj Remark 10.4]. These stabilization pa- 
rameters may depend possibly also on time steps for time-dependent problems, so their choices 
become more tricky in practice as we have to balance between temporal and spatial errors [8]. 

By changing the sign of the convective term in the weighted least-squares formulation, the 
unusual stabilized finite element method (USFEM) can achieve the absolute stability for any 
positive stabilization parameter involved in the scheme, but it is still a tricky and inconclusive 
technical issue of how to choose this parameter in order to obtain good accuracy [9l 110^ [TT | I12j. 
The variational multiscale method was developed based on the inherent multiscale structure of 
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solutions [121 HH ISl in] • This method defines the large scales by a projection into an appropriate 
subspace, but may also involve the technical issue of how to select a stabilization parameter to 
balance the stability and accuracy. 

As it is known [5], explicit Galerkin solutions for flow problems could be quite under-diffusive, 
effectively increasing the Peclet or Reynolds number. Furthermore, explicit methods are gener- 
ally conditionally stable. But explicit schemes have the advantages that they may not need to 
solve systems of algebraic equations |17j or the resulting stiffness matrices stay the same in the 
time marching process. 

The characteristic-based-split (CBS) method has been widely studied for fluid and solid 
dynamic problems \18\ \T9\ [201 [2T| , and we refer to the monograph [22] and the references therein 
for its detailed introduction and various applications. This method is based on the splitting of 
the convection and diffusion parts. The convection part is formally handled by the standard 
characteristics method, where the numerical solutions at the current time are updated by the 
approximations at the previous time. But the schemes need to locate some spatial points based 
on the characteristics, and the spatial points are likely no longer grid points of the spatial 
discretization. One way to avoid this is to adjust the meshes, while another way is to apply the 
standard interpolation to evaluate the values of the solutions at these spatial points using the 
values of the solutions and other quantities at grid points. An alternative technique, used in the 
CBS method, is to approximate numerical solutions at computed spatial points by the solutions 
and other quantities at grid points by Taylor expansion. In addition, the CBS method needs to 
approximate the average convective field, for which different treatments may lead to different 
schemes, such as fully explicit, semi- implicit or implicit ones, and also different stabilization 
effects [211 [22]. 

In the derivation of the new scheme, we shall use the same operator splitting as the CBS 
method did, to split the convection diffusion system into a purely convective part and a diffusive 
part. The diffusion part is discretized by the standard backward scheme. But the central 
difference from the CBS method lies in our new treatment of the convection part, which is 
completely independent of the characteristic curves and any spatial grid points used, unlike the 
CBS method. 

Another novel idea of the new method is the flexibility in its special explicit treatment of the 
convection part: we can recursively execute the explicit convection step up to a finite number of 
times with smaller local time steps during one diffusion correction. This can essentially improve 
the stability of the resulting scheme so it may behave like an unconditionally stable scheme. 

The rest of the paper is arranged as follows. The single-step scheme is first derived for the 
convection diffusion equation in Section [27T1 and its multistep variant in Section [2r2l The new 
scheme is then extended in Section[3| for the Navier-Stokes equations. Numerical experiments 
are carried out in Section[Hto check the accuracy, stability and performance of the new schemes, 
as well as to investigate how the stability condition can be improved by the multistep scheme 
compared with the single-step one. At the end of this numerical section, the driven cavity flow 
problem is tested with the new scheme and compared with the benchmark results to demonstrate 
the validity of the new method. Some concluding remarks are given in Section[5l 
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2 Derivation of algorithms 

In this section we derive a new method for solving the convection-dominated diffusion equation 
([1]). For the purpose we introduce some notations. We first partition the time interval [0,r]: 
= to < ti < ■■■ < tN = T, with tn = nAt and At = T/N. We will use u"- and u'^^^ 
respectively for the approximate values of u{-,t) at t = t„ and t„ + At/2. But when u{-,t) is 
a known function, and n"'^2 will stand for its exact values at t = t„ and tn + At/2, e.g., 
r = /(•,*„), and b« = b(-,t„). 



2.1 Single-step scheme for the convection diffusion equation 

We first adopt the standard operator splitting technique [3] and split the convection diffusion 
equation ([1]) into a pure convection equation and a diffusion equation. Then we approximate two 
equations in time by the central difference and backward Euler schemes respectively to obtain 



At 



+V-(b"+2n"+2) = (5) 
- V • (eVn"+i) + c'^+in'^+i = g''+\ (6) 



where / and g can be any functions such that F = f + g. However in order to have a unified 
principle for the selection of the components / and g for both the convection diffusion equa- 
tion and Navier-Stokes equations, we will suggest some special choice of / and g later on; see 
Remark l3.1l 

We shall use finite element methods to solve ([5]) and ([6]) respectively for the solutions u^'^^ 
and tt"+^ To do so, we need the variational formulations of these two equations. For equation 
([6]), it is straightforward to derive its variational form: 

Find u"+i € H^{Vt) such that = on F and solves 



"+\7;) + At(eVtx"+Sv?;) + At(c"+^n"+Si;) = «+\ ?;) + At(5"+\ ^ v e hI{^1) . (7) 



On the other hand, the solution of the convection step ([5]) is more tricky. Clearly the scheme 
is implicit and involves the solution of a linear convection equation. The main idea of this work 
is to propose an explicit scheme to solve this linear convection equation. For this aim, we apply 
the Taylor's expansion to compute n"^2 by the values at previous times, and can write 



u'^+s ^ n(x, tn + ^) = u(x, tn) + ^nt(x, tn) + 0{At^), 
then using the convection equation 

ut + V-{hu) = f (8) 

we deduce 

u"+3 + V-(bV)) =:^". (9) 

Using this relation, we can rewrite ^ as 

,,n+l _ / , 1 \ ,1 

+ V- (b"+2^") = /"+2. (10) 



At 
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Noting that ([S]) is a pure convective equation, only partial boundary condition on the inflow 
boundary should be imposed, namely 

:= {x G T; b(x, t) ■ n(x) < 0} (11) 

where n(x) is the outward normal to the boundary of at x. Accordingly we should set a 
similar condition on the inflow boundary associated with the scheme (jlOp . So for any positive 
integer n, we define 

r- :={xGr; b"(x)-n(x) <0}. (12) 

As the exact solution is specified on the entire boundary (cf. ([I])), it is natural for us to assume 
the values for the solution u^~^^ to pU|) on the inflow boundary 

on r-^,. (13) 



This induces the following test space for the scheme (jlOp : 

i^i (O) = {wG H\n)- w = o on r;+i} . 

Now multiplying a test function v £ (J7) on both sides of (|lUp . and integrating over Q and 

-'- n + l 

using the integration by parts we obtain 

{u:+\v) = {u^,v) + At{f^-^-2,v) 

+ At{C, b"+^ . Vv) -At<C, t-b^+s . n >r\r-+, 
= («",?;) + At(/ "+5,?;) 

+ At (^n" + ^ (/" - V • (b"u")) , b"+5 . Vv^ (14) 

- At(u" + — (/"- V • (b"n")) , 7;b"+^ • n)p> p- . 

It remains to introduce the spatial discretizations for both equations d?]) and (|14p , which we 
will do by finite element methods. Assume that Vh is a finite element space approximating the 
Sobolev space H^{Q,), and Ih is the interpolation operator of H^{^) into Vh- Then based on the 
variational formulations (I14p and ([7]), we propose the following single-step scheme for solving 
the convection-dominated diffusion problem ([1]). 

Algorithm 1 (Single-step scheme). 
Step 0. Compute the initial value = IhUQ. For each n = 0, 1, • • • , — 1, do the following. 
Step 1. Find G Vh such that m^+^ = IhUb~^^ on r~_^_^ and it solves 

Kf,Vh) = {ulvh) + Atif-+Kvh) 



At 



At(u^^ + ^(r-V- (b"nD ) , b"+^ . v^. 
At 



At(ul + (/- - V • (b"^x^)) , Vhh^+'^ • n) yvhC^Vhn H^^ {Q) . 

\ 2 V / /r\r in+i 
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Step 2. Find € Vh such that = 4^^"^^ on T and it solves 

Remark 2.1. One may compute the term {u'^']^^,Vh) in Step 1 by the standard mass-lumping 
technique JT^ , then u^^^ can be computed explicitly without solving a linear system. 

2.2 Multistep scheme for the convection diffusion equation 

The focus of this work is mainly on the case when the convection diffusion system ([T|) is 
convection-dominated. For this case the stability of the explicit single-step scheme (Algorithm 
1) may pose severe restrictions on time steps, leading to sufficiently small time steps and great 
computational efforts for the entire numerical resolution process. 

To improve the stability, we may execute the convection step (Step 1) a few times for each 
diffusion correction (Step 2) so that we can use much smaller time steps for the convection part 
and much larger time steps for the diffusion step. To do so, we write the result u^'\^^ of Step 1 
formally as 

<f = F^,z{^t,r,r+\h-,h-+\uiu^+') . (15) 

Then the multistep scheme is to run this convection step m times with smaller time step size 
At/m for , namely we compute 

= ^^(^>r"^>r"->b-^>b"^->C^>-r-), (16) 

recursively for i = 1, 2, • • • ,m, with u]^^ = u^- 

We shall call 6t = At/m and At as the local time step size and the global time step size 
respectively. Replacing Step 1 by the multistep iteration (jl6p , we propose the following multistep 
scheme for the convection diffusion equation ([T]). 

Algorithm 2 (Multistep scheme vi^ith index m). 
Step 0. Compute the initial value = IhUq. For each n = 0, 1, • • • , — 1, do the following. 

n+— n-\-— n+ — 

Step 1. Set ^ = u^. For i = 1,2, ■ ■ ■ ,m, compute % ^"^ ^ Vh such that u^^™- = IhU^^ on 
^n+j/m solves for all Vh GVhD H^^ {^), 

n-\-i/m, 

- St (ui::^ + I (/-^ - V . (b-^n--)) ,.,b-^ . n)^ ^ . 

"^"^"^ri + i/ m 

Step 2. Compute n^+^ G Vh such that u^^^ = hu^^^ on P and it solves for ah Vh G V/inFo(r2), 
iul+\vh) + At{eVul+\Vvh) + At{c^+\l+\vh) = {ul+\vh) + At{g^+\vh) . 
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3 Single-step and multistep schemes for Navier-Stokes equa- 
tions 

We are now going to extend the new schemes proposed in Sections l2. 1112.21 for the convection- 
dominated diffusion equation to the Navier-Stokes equations ([3]). For the purpose, we spht the 
system ^ into a pure convection system and a diffusion system (the generalized Stokes problem) 
as fohows: 



At 



■ V)u"+^ 


= f"+3 


+ Vp"+^ 


= g"+^ 


V • u"+i 


= 0. 



(17) 



, n+l _ 

i?e-iAu"+^ + Vp"+^ = g"+\ (18) 

(19) 

It is straightforward to derive the variational form of the generahzed Stokes system (jl8p - (jl9p : 
Find u''+^ G H^(0) and p G Lg(0) such that u'^+i = u^+^ on F and it solves 

(Af)-i(u"+i,v)+/2e-i(Vu"+SVv)-(p"+\V-v) = {AtrHu",+\^r) + {g^-^\^r) , (20) 

(V-u"+Sg) = (21) 

for any v G Ho(r2) and q G LI{Q). 

Next we will do the same as we did in Section [2?T] to propose an explicit scheme for solving 
the convection system (jl7p . To do so, we first handle the nonlinear convection term involving 
In fact, combining the Taylor's expansion 

u"+5 ^ u(x, tn + ^)= U(x, tn) + ^M^, ^n) + O(Ai^), 

and the pure convection equation 

ut + {u- V)u = f , (22) 
we can obtain a similar approximation to ^ but in a vector-valued form: 

u"+5 ~ u" + ^ (r - (u" • V)u") =: r?". (23) 

Again, we introduce the inflow boundary 

r-+, = {x G fl; u^+i • n(x) < O} . 
Then we can write by using integration by parts for any v G H^(il) with vL- = that 

((?f • V)?f , v) = (r?", 7f . n v)„ p- - (7f , V • 7?"v) - (r?", (r?" • V)v) , (24) 

using this relation and plugging (j23|) in (fT7|) we derive the variational form of (pT|) : 

«+\v) = (u",v) + At(r+iv)-At(r?",(77".n)v) 

+Ai(r?",(V-ry" +77"-V)v) . (25) 
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Remark 3.1. We observe from the formulation i25\) that V • f " is needed in the term V • r/", 
hence it adds some extra regularity on the source component f . This suggests us to better choose 
f = m the decomposition F = f + g for the Navier-Stokes equations so that the new scheme 
does not need the evaluation o/V-f", unlike in the CBS methods \2(^ \21\ [2 ^ . For the unification 
of the numerical schemes for both the convection diffusion equation and Navier-Stokes equations, 
we shall always select f = in Algorithms 1 and 2 from now on for the convection diffusion 
equation. 

Let V/j and be two finite element spaces approximating the Sobolev space H^(r2) and 
Lq{Q), and Ih be a interpolation operator of H^(il) into V/j. By virtue of the variational formu- 
lations (I25p and (I20p . we propose a single-step scheme for solving the Navier-Stokes equations 
©• 

Algorithm 3 (Single-step scheme). 
Step 0. Compute the initial value = I/iUq. For each n = 0, 1, • • • , — 1, do the following. 
Step 1. Find u^+^ G Vh such that u)J+^ = IftU^"^^ on r~_^^ and it solves 

K,v,)-A^(<,(r?^nK) 

^ ' J- „+i 

+At{r^Ji, {V-r,l + th- V)v^) Vv„ G V;, n ■ 

Step 2. Find u^+^ G and ph e Mh, such that u^+^ = Ihu"^-^ on F and it solves 

At-^K+i, V,) + Re-HVul+\Vv,) - {pl^\V ■ v,) = At-Hul+\v,) + (g"+\ v,), 

iV-nl+\qf,) = 

for any v/, G n Hj(J]) and qn ^ Mh- 

For simplicity we have used in Algorithm 3 the notation r/^, which is defined as in (|23p 
but with u" replaced by u^. Similarly we shall use the following notation in Algorithm 4: 

We observe from Algorithm 3 that the nonlinear convection term (u • V)u in Navier-Stokes 
equations has been treated explicitly in the time marching process, which may severely restrict 
the time step size in order to ensure the stability of the scheme when the convection is dominated 
in comparison with the diffusion of the flow system. To improve the stability, we may apply 
Step 1 several times with a smaller time step size during one diffusion correction (Step 2). For 
this purpose we write the result of Step 1 formally as 

Kt.' = Ffi, (At, (26) 

Then a multistep variant of this scheme is to execute this step m times with a smaller time step 
size At/m to derive u^^^: 

^h,r = Ffi„(— ,u^,^'",u, (27) 



K:t^v) = 
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for i = 1,2,- ■■ ,m, with u^^ = u^. This leads to the fohowing multistep scheme for the 
Navier-Stokes equations. 

Algorithm 4 (Multistep scheme with index m). 
Step 0. Compute the initial value = I/iUq. For each n = 0, 1, • • • , — 1, do the following. 
Step 1. Set ^ = u^; then for i = 1, 2, • • • ,m, 

n+— n+— n+— _ 

compute ^™ G \ ^ such that ^"^ = i/jUj^ " on i and it solves 

/ \ / n+— . / n+— , n+— > \ 

\ /r\r- ,, 

+'^ifC'^'(^-C^ + C^-^)^0 Vv,GV,nHi_ (0). 

\ J n-\-ijm 

Step 2. Compute (u^+Sp^+^) G V/, x M/, such that u^+^ = I/iU^+^ on V and it solves 
At-^K+i, v;,) + i?e-i(Vu;:+\ Vv/,) - V • V,) = At-i(u;:+\ v;,) + (g"+\ V,), 

(v.<+\g.) = 

for any (v/,, g/,) G (V/, n HJ(J7)) x M/,. 

Remark 3.2. The second steps in Algorithms 3 and 4 can be replaced by the projection-type 
methods so that the pair of finite element spaces for approximating the velocity and pressure does 
not need to meet the LBB condition and only Poisson problems are needed to solve for updating 
both the velocity and pressure. For the projection method, we refer to the pioneering work by 
Chorin JMI and Temam 



4 Numerical experiments 

In this section we shall carry out two sets of numerical tests to check the actual convergence 
orders of the single-step and multistep schemes proposed in the previous two sections and how 
the multistep scheme improves the stability of the single-step scheme. 

Let 7h be a regular triangulation of 17, with hx = diam(i^) for K G 7^, and h = max^'eTj^ hx- 
We shall use the following linear finite element space Vh C H^{Q): 

Vh = {wh £ H\ny, Wh\K e Pi{K) ^KeTh} (28) 

for the solution of the convection diffusion equation (JT]), and the following Taylor- Hood finite 
element spaces [25] 

V;, = {^rh£H\nf;vh\KC^P2{Kf yKGTh}, (29) 
Mh = {qh e H\ny, qh\K e Pi{K) yKeTh} (30) 

for the solution of the Navier-Stokes Equations ([3]). 

We recall that we have used the central finite difference scheme for the convection diffu- 
sion equation and the backward Euler scheme for the diffusion equation in time discretization. 
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Therefore it is natural for us to expect the following numerical convergence orders when the 
finite element spaces in ([28]) and (j30|) are used: 

Wu"" -u^h2(^n)<C{h^ + At) 

for the convection diffusion equation ([T|), and 

I|u^-u;^||i2(^) <C(/i3 + At) and Wp"" - ^2^^) < C{h^ + At) 

respectively for the velocity and pressure of the Navier-Stokes equations ([3]). Naturally we may 
think that the convergence of the scheme can be improved by using a second-order scheme 
for the diffusion equation in time discretization, but our numerical experiments have firmly 
disapproved this conjecture. We are currently investigating the possible treatments to help 
construct the schemes which have second order temporal accuracy. 

We remark that all the errors shown in this section are the L^-norm errors at the terminal 
time t = T unless specified otherwise. 

4.1 Tests for the convection diffusion equation 

We first apply the new single-step and multistep schemes to the following two examples which 
are taken from references [8] and |16] . 

Example 1. The coefficients and domain in equation ([2]j are taken to be the following: 

d = 2, T = l, e = 10~^ b = (l,-l)^, c=l, n = {0,lf 

with the exact solution given by u{x,y,t) = e^'^* sin(27rj;) sin(27ry). 

This example is a slight modification of the one in [8], where e*^™^^'^*) is used. Instead we use 
e^'^*, which makes the solution vary in a much larger range, namely in the interval [— e^'^, e^'^], 
and has a much larger norm, \\u{-, 1)|| = ^e^'^ ~ 267.7458. 

Example 2. The coefficients and domain in equation (Op are taken to be the following: 

d = 2, T = l, e = 10"^ b = (2, -1)^, c = l, n = {0,lf 
with the exact solution given by u{x,y,t) = t^ cos{xy'^) . 

To compute the actual convergence orders of the numerical schemes, we shall use the uniform 
triangulations of domain Vt with triangular elements in all our numerical simulations. 

4.1.1 Convergence Tests for the single-step scheme 

In order to find the actual convergence order of the single-step scheme (Algorithm 1) in time, we 
choose a very small mesh size and then observe the changes of the errors when the time step size 
is halved. Similarly we will do the other way around when we try to find the actual convergence 
order of the single-step scheme (Algorithm 1) in space. 

Tables 1 and 2 show the L^-norm errors with different mesh sizes when the time step size is 
fixed for Examples 1 and 2 respectively. Clearly we see the second order spatial convergence of 
the single-step scheme (Algorithm 1). 
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Table 1: Convergence results of Algorithm 1 for Example [T] with fixed At =1/2 



h \\u — n/j|| order 

1/4 9.18526(+1) 

1/8 1.84780(+1) 2.3135 

1/16 4.24054 2.1235 

1/32 1.03466 2.0351 

1/64 2.54797(-l) 2.0217 

1/128 6.36669(-2) 2.0007 



Table 2: Convergence results of Algorithm 1 for Example [2] with fixed At = 1/2 



h \\u — Uh\\ order 

T/l 9.76826(-3) 

1/8 2.41756(-3) 2.0145 

1/16 6.02478(-4) 2.0046 

1/32 1.49729(-4) 2.0086 

1/64 3.69132(-5) 2.0201 

1/128 8.70186(-6) 1.9687 



Now we fix the uniform mesh size at /i = 1/128, and run the single-step scheme (Algorithm 
1) for Examples [T] and [2] with the following sequence of time step sizes 

At = 0.1/2'=, /c = -1,0,1, 2,--- (31) 

to find out the stability region of the numerical scheme. The numerical results are listed in Tables 
[3]and[31 from which we observe that Algorithm 1 does not converge till k = 6 and 7 respectively 
for Examples [1] and [21 corresponding to two rather small time step sizes of At = 1 /640 and 
1/1280. Such restrictions on time step sizes are natural, required by the stability condition for 
the explicit time marching scheme we have used for the convection step in Algorithm 1. As we 
shall see in the next subsection, the new multistep scheme can essentially improve the stability 
condition. 



Table 3: Convergence results of Algorithm 1 for Example [T] with fixed h = 1/128 



At 


\\u - UhW 


order 


0.1/2^ 


divergence 




0.1/2^ 


2.02151 




0.1/2'^ 


1.01181 


0.9985 
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Table 4: Convergence results of Algorithm 1 for Example [2] with fixed h = 1/128 



At 



U-Uh 



order 



0.1/2' 
0.1/2' 



divergence 
1.71484(-4) 



,7 



4.1.2 Stability improvement by the multistep scheme 

We can observe from the previous subsection that the single-step scheme (Algorithm 1) may pro- 
vide the expected convergence and preserve the accurate convergence orders when it converges. 
However, this scheme requires sufficiently small time step size as shown in Tables [3] and [H hence 
may restrict its applications in practice. The multistep scheme (Algorithm 2) is proposed to 
improve the stability of the single-step scheme. This section is to test how the multistep scheme 
can improve the stability region. 

We note that At is the global time step size, which is used for the diffusion correction. As 
we are interested mainly in the convection-dominated diffusion problems, the time step size 
required for the convection is usually much smaller than the one for the diffusion. 

In our numerical tests, for each fixed global time step At = 0.1/2^^ (k = —1, 0, 1, 2, • • • ), we 
run the multistep scheme with index m = 1,2^,2^,-- - until we observe the convergence of the 
scheme, and then record the corresponding index m; see Tables [5] and [6] for the recorded index 
m corresponding to each fixed At and the resulting relative L^-norm error of the approximate 
solution. 

As we see from Table [5l when we take At = 0.1, which is too large for the stability of 
the explicit scheme involved in the convection step, but we can still achieve the convergence 
of the multistep scheme with index m > 64. Tables [5] and [6] have demonstrated that though 
the single-step scheme does not converge for a fixed global time step At, the multistep scheme 
always converges when the index m is appropriately large. So we can conclude that if we 
take an appropriately large index m, say m = 30, the multistep scheme can be viewed as an 
unconditionally stable scheme. 

Furthermore, we have also computed the convergence orders of the multistep scheme in 
terms of the global time step size for Examples 1 and 2 with a fixed index m and mesh size 
h. The results are shown in Tables [7] and [8l Combining these results with the ones for the 
single-step scheme (cf. Table [3]), we can clearly observe the first order temporal convergence for 
both examples. 

Next, we carry out some numerical tests to check how the multistep scheme can improve the 
stability region quantitatively. For each fixed mesh size h, we increase the index m gradually 
and record the largest global time step size At that can ensure the convergence of the entire 
algorithm. And the largest time step size will be written as the critical time step size Atcrit 
for the stability of the algorithm. The results are shown in Tables [9] and [TOl from which we 
can see that the stability region is nearly doubled when the index m of the multistep scheme is 
doubled. So the multistep scheme can indeed clearly and essentially enlarge the stability of the 
entire algorithm. 

We remark that we have done many more numerical experiments for Examples 1 and 2, 
but with the diffusion coefficients e varying in a wider range, from 

10-3 to 10-l^ and many 
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Table 5: 


Stability of Algorithm 2 for Example [T] with index m and fixed h 


= 1/128 


At 


m 


\\u\\ 


0.1/2'^ 

0.1/2^ 

0.1/2^ 

0.1/23 

0.1/22 

0.1/2^ 

0.1 

0.2 


1 

2 

4 

8 
16 
32 
64 
128 


7.55011(-3) 
1.54379(-2) 
3.15569(-2) 
6.41671(-2) 
1 30693f-l') 
2.69788(-l) 
5.68483(-l) 
1.26837 


Table 6: 


Stability of Algorithm 2 for Example [2] with index m and fixed h 


= 1/128 


At 


m 


\\u\\ 


0.1/2'' 

0.1/2^ 

0.1/2^ 

0.1/2^ 

0.1/2^ 

0.1/2^ 

0.1/2^ 

0.1 

0.2 


1 
2 

4 
8 

16 
32 
64 
128 
256 


1.65924(-4) 
7.46950(-4) 
1.96066(-3) 
4.39337(-3) 
9.28328(-3) 
1 96240f-2) 
4.11452(-2) 
8.41251(-2) 
1.70957(-1) 


Table 7: Convergence order of Algorithm 2 for Example [T] with fixed index m = 64 and h = 


At 


11^ - UhW 


order 


0.1 

0.1/2^ 
0.1/22 
0.1/2^ 
0.1/2^ 
0.1/2^ 
0.1/2^ 


1.52209(+2) 

7.23099(+l) 

3.51127(+1) 

1.73718(+1) 

8.66679 

4.34772 

2.18766 


1.0738 
1.0422 
1.0152 
1.0032 
0.9952 
0.9909 
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Table 8: Convergence order of Algorithm 2 for Example [2] with fixed index m = 128 and h = 



At 






ll'" - UhW 




order 


0.1 






8.69440(-2) 






0.1/21 






4.27948(-2) 




1.0227 


0.1/2^ 






2.07747(-2) 




1.0426 


0.1/2^ 






1.00067(-2) 




1.0538 


0.1/2^ 






5.01775(-3) 




0.9959 


0.1/2^ 






2.54242(-3) 




0.9808 


0.1/2^ 






1.31725(-3) 




0.9487 


Table 9: Critical ^ 


rlobal time step 


size At, 


crit of Algorithm 2 for Example [T] 


in terms of index m 


m 


1 


2 


10 


20 


40 80 


h = 1/64 












Atcrit 


0.0049 


0.0093 


0.046 


0.093 


0.18 0.37 


h = 1/128 












^^crit 


0.0024 


0.0045 


0.022 


0.045 


0.091 0.18 


Table 10: Critical 


global time step size Atcrit of Algorithm 2 for Example [2] 


in terms of index m 


m 


1 


2 


10 


20 


40 80 


h = 1/64 












Atcrit 


0.0032 


0.0060 


0.029 


0.058 


0.11 0.23 


h = 1/128 














0.0015 


0.0030 


0.014 


0.028 


0.057 0.11 
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different convective vectors b, and observed similar convergence and stability behaviors for the 
single-step and multistep schemes as we have shown above. 

4.2 Tests for the Navier-Stokes Equations 

Now we will apply our new single-step and multistep schemes (Algorithms 3 and 4) to two 
examples of Navier-Stokes equations with analytical solutions to check the actual convergence 
orders of the schemes and how the multistep scheme improves the stability of the single-step 
scheme. Then we will apply these schemes to the benchmark problem of the lid-driven cavity 
flow to verify their validity. 

Example 3. Consider the Navier-Stokes equations ^ with the following parameters: 

O = [0, l]^ T=l, Re = 5000 and 10000 

with the exact solution (u,p) = {ui,U2,p) given by p = (x^ — y^)cos(t) and 

ui = 10x^{x - lfy{y - l){2y - 1) cos(t) , U2 = -I0x{x - l)(2x - l)y'^{y - if cos(t) . 

Example 4. Consider the Navier-Stokes equations ^ with the same parameters as in Exam- 
ple\^ hut the exact solution (u,p) = {ui,U2,p) given by 

ui = t^y"^ , U2 = t^x, p = tx + y - {t+l)/2. 

This is an example where only a discretization error in time occurs 126}/ . 

4.2.1 Convergence Tests for the single-step scheme 

We first verify the convergence orders of the single-step scheme (Algorithm 3) in both space 
and time for Example [3j Tables [TT]fT2] present the convergence results in time for the Reynolds 
numbers Re = 5000 and 10000 respectively, with a fixed uniform mesh of size h = 1/128, and 
Tables [T3lfHl give the convergence results in space for the Reynolds numbers Re = 5000 and 
10000 respectively, with a fixed At = 10~^. From these tables we can clearly see the optimal 
first order convergence of the single-step scheme in time and the optimal third and second order 
convergence in space respectively for the velocity and pressure. 

For Example IU we have tested the single-step scheme (Algorithm 3) with the Reynolds 
numbers Re = 5000 and 10000, and the uniform mesh of size h = 1/48 and 1/64, and the 
sequence of time step sizes as listed in (131 p . The results have shown that the scheme converges 
only when the time step size At = 0.1/2'^ is sufficiently small, namely when k takes at least 4 
{At = 1/160) and 5 {At = 1/320) respectively for /i = 1/48 and 1/64. This test indicates that 
the single-step scheme may require sufficiently small time step size to ensure its convergence, as 
one can expect for this strongly convection-dominated example. In the next Section l4.2.2l we will 
show the multistep scheme (Algorithm 4) can essentially improve the stability of the single-step 
scheme. 

4.2.2 Stability improvement by the multistep scheme 

As shown in the previous subsection, the convergence of the single-step scheme (Algorithm 3) 
for Example m requires a sufficiently small global time step size for a fixed mesh size h. In order 
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Table 11: Convergence of Algorithm 3 for Example [3] with h = 1/128 and Re = 5000 



At 




order 


\\P-Ph\\ 


order 


0.2 


3.28203(-3) 




1.00222(-4) 




0.1 


1.65607(-3) 


0.9868 


4.79084(-5) 


1.0648 


0.1/21 


8.31889(-4) 


0.9933 


2.35900(-5) 


1.0221 


0.1/22 


4.16919(-4) 


0.9966 


1.20411(-5) 


0.9702 


Table 12: 


Convergence of Algorithm 3 for Example [3] with h = 1/128 and Re = 


10000 


At 


l|u - ^h\\ 


order 


\\P-Ph\\ 


order 


0.2 


3.28203(-3) 




9.98106(-5) 




0.1 


1.65607(-3) 


0.9872 


4.77006(-5) 


1.0652 


0.1/2^ 


8.31889(-4) 


0.9935 


2.34866(-5) 


1.0222 


0.1/2^ 


4.16919(-4) 


0.9967 


1.19910(-5) 


0.9699 



Table 13: Convergence of Algorithm 3 for Example [3] with At = 10"^, Re = 5000 and T = 0.2 



h ||u — u/i|| order order 

1/4 1.31468(-3) - 6.45683(-3) 

1/8 1.81020(-4) 2.8607 1.61419(-3) 2.000016 
1/16 2.38018(-5) 2.9270 4.03547(-3) 2.000002 
1/32 3.00134(-6) 2.9874 1.00887(-4) 1.999996 
1/48 8.71038(-7) 3.0511 4.48386(-5) 2.000004 



Table 14: Convergence of Algorithm 3 for Example [3] with At = 10"^, Re = 10000 and T = 0.2 



h 


l|u - Uh\\ 


order 


\\P-Ph\\ 


order 


1/4 


1.31577(-3) 




6.45686(-3) 




1/8 


1.81589(-4) 


2.8572 


1.61419(-3) 


2.000016 


1/16 


2.42194(-5) 


2.9064 


4.03547(-4) 


2.000002 


1/32 


3.20061(-6) 


2.9197 


1.00887(-3) 


1.999996 


1/48 


9.28064(-7) 


3.0533 


4.48386(-5) 


2.000004 
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to improve this severe restriction on time step size by the single-step scheme, we now show how 
we can achieve the convergence for large global time step size by the multistep scheme. For each 
fixed At = 0.1/2^^' (k = —1, 0, 1, 2, • • • ), we run the multistep scheme with index m = 1, 2^, 2^, • • • 
until we observe the convergence of the scheme, and then record the corresponding index m; 
see Tables [15] and [16] for the recorded index m corresponding to each fixed At and the resulting 
relative L^-norm errors of the approximate solutions for the velocity and pressure. 

As we see from Table [T5] when we take the global time step At = 0.1, which is too large for 
the stability of the explicit scheme involved in the convection step, but we can still achieve the 
convergence of the multistep scheme with index m > 32. Tables [T5] and [TBI have demonstrated 
that though the single-step scheme does not converge for a fixed At, the multistep scheme 
always converges when the index m is appropriately large. So we can conclude that if we 
take an appropriately large index m, say m = 30, the multistep scheme can be viewed as an 
unconditionally stable scheme. 

Next we have tested the actual convergence orders of the multistep scheme when the index 
m is fixed at m = 64. Tables [TTlfTS] have showed the computational results for Re = 5000 and 
10000 with fixed h = 1/48 and 1/64 respectively. We can observe clearly the optimal first order 
convergence for both velocity and pressure in terms of the global time step size. 

The last test we have carried out is to check how the multistep scheme can improve the 
stability region quantitatively. For each fixed mesh size h, we increase the index m gradually 
and record the largest global time step size At (the critical time step size Atcrit as we called 
earlier) that can ensure the convergence of the entire algorithm. The results are shown in Table 
[T9] from which we can see that the stability region is nearly doubled when the index m of the 
multistep scheme is doubled. So the multistep scheme can indeed clearly and essentially enlarge 
the stability of the entire algorithm. 

We end this subsection with some concluding remarks on convergence and stability behaviors 
of the single-step and multistep schemes, based on our observations from the numerical tests in 
this and previous subsections. 

• The single-step scheme (Algorithm 3) is generally conditionally stable, and requires suffi- 
ciently small time step size to ensure its convergence with a fixed mesh and larger Reynolds 
number. 

• The multistep scheme (Algorithm 4) can essentially relax the restriction of the time step 
size (see Tables [T5l 1161 and I19p . and behaves like a nearly unconditionally stable scheme. 

• Comparing the results in Tables [T5lfT6] with the ones in Tables [T71ll8|, we can clearly see 
the stability and robustness of the multistep scheme (Algorithm 4). For example, for the 
global time step size At = 0.1/2^, the multistep scheme with a small index like m = 2 and 
a large index like m = 64 provides about the same accuracies; see Tables [16] and [THl 

4.2.3 The lid-driven cavity flow 

As our final numerical example we test a popular benchmark problem, i.e., the lid-driven cavity 
fiow problem, where the fluid is enclosed in a unit square box, with an imposed velocity of unity 
in the horizontal direction on the top boundary, and a no-slip condition on the remaining walls. 
We shall compare our results with three benchmark results: Ghia et al. [27] with h = 1/128 for 
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Table 15: Stability of Algorithm 4 for Example U] with index m and fixed h = 1/48, Re = 5000 

At m ||u - Ufell \\p-Ph\\ 

0.1/24 I 1.61352(-3) 9.22561(-3) 

0.1/23 4 3.15065(-3) 1.82083(-2) 

0.1/22 8 6.32378(-3) 3.52747(-2) 

0.1/2^ 16 1.34098(-2) 6.64021(-2) 

0.1 32 2.70632(-2) 1.18921(-1) 

0.2 64 5.60465(-2) 1.97385(-1) 



Table 16: Stability of Algorithm 4 for Example H] with index m and fixed h = 1/64, Re = 10000 



At m ||u - UhW \\P-Ph\\ 

0.1/2^ 1 8.13177(-4) 4.65014(-3) 

0.1/2^ 2 1.61134(-3) 9.24510(-3) 

0.1/23 4 3.46105(-3) 1.82018(-2) 

0.1/2^ 16 6.37609(-3) 3.52913(-2) 

0.1/2^ 32 1.29496(-2) 6.64033(-2) 

0.1 64 2.67569(-2) 1.18926(-1) 

0.2 128 5.67472(-2) 1.97454(-1) 



Table 17: Convergence order of Algorithm 4 for Example H] with fixed h = 1/48, Re = 5000 and 
fixed index m = 64 



At 




order 


\\P-Ph\\ 


order 


0.2 


5.60465(-2) 




1.97385(-1) 




0.1 


2.64298(-2) 


1.0845 


1.18888(-1) 


0.7314 


0.1/2^ 


1.28200(-2) 


1.0438 


6.63946(-2) 


0.8405 


0.1/22 


6.33073(-3) 


1.0180 


3.52954(-2) 


0.9116 


0.1/2^ 


3.17694(-3) 


0.9947 


1.82354(-2) 


0.9527 


0.1/2^ 


1.60997(-3) 


0.9806 


9.27332(-3) 


0.9756 


0.1/2^ 


8.20838(-4) 


0.9719 


4.67745(-3) 


0.9874 
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Table 18: Convergence order of Algorithm 4 for Example H] with fixed h = 1/64, Re = 10000 
and fixed index m = 64 



At 




order 




\\P-Ph\\ 




order 


0.1 


2.67569(-2) 






1.18926(-1) 






0.1/21 


1.29252(-2) 


1.0497 




6.64054(-2) 




0.8407 


0.1/2^ 


6.35955(-3) 


1.0232 




3.52978(-2) 




0.9117 


0.1/2^ 


3.17773(-3) 


1.0009 




1.82330(-2) 




0.9530 


0.1/2^ 


1.60477(-3) 


0.9856 




9.27095(-3) 




0.9758 


0.1/2^ 


8.16528(-4) 


0.9748 




4.67451(-3) 




0.9879 


Table 19: Critical global time step size 


Atcrit of Algorithm 4 for Example U] in 


terms of index m 


m 


1 


5 


10 


20 


40 


80 


Re = 10000, h 


= 1/64 














0.0039 


0.018 


0.024 


0.048 


0.089 


0.18 



Reynolds numbers Re = 100,400, 1000 and 3200; Erturk et al. [28] with h = 1/128 for Reynolds 
number Re = 1000; Botella et al. [22] for the Reynolds number Re = 1000. 

In all our computations for this example, we use the uniform mesh of size h = 1/128 
and the Taylor-Hood elements ()30p . and have tested the cases with Reynolds numbers Re = 
100, 400, 1000 and 3200, and the global time step size At = 0.004. The stoping condition for 
time advancing, which is considered as the criterion of capturing the steady state solution, is 
chosen as 



\^h 



n+1 1 



< 10"^ 



where is the finite element solution at time t = tn- We have observed from our numerical 
results that the single-step scheme (Algorithm 3) works when the Reynolds number is relatively 
small, e.g.. Re = 100,400 and 1000, but it is unstable when Re is large, e.g.. Re = 3200. But 
the multistep scheme may still work for larger Reynolds number, e.g.. Re > 3200. 

Tables [20]|21l present the streamfunction values and the locations of the primary and sec- 
ondary vortices for various Reynolds numbers. Figures [H [2] and [3] show the computed velocity 
components and vorticity profiles along the horizonal and vertical lines compared with the re- 
sults of Ghia et al. [27] and Botella et al. [29]. As one can see that the results by the new schemes 
confirm very well with the ones by the benchmark schemes. 



5 Concluding remarks 

We have proposed a new splitting method for solving time-dependent convection-dominated dif- 
fusion problems. A pure convection problem and a pure diffusion problem are solved successively 
at each iteration of the method. Explicit schemes are proposed for the time discretization of 
the convective problem. The explicitness of the scheme may cause a severe restriction on the 
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Table 20: Streamfunction values ^min, ^max and locations of the primary and secondary vortices 



Vortex 


property 




Re=1000 


Re=1000 


Re=1000 








Single-step scheme 


Ghia et al. [27] 


Erturk et al. [28] 


Primary 






-0.114722 


-0.117929 


-0.118781 




Location (x, 


y) 


(0.5313, 0.5625) 


(0.5313, 0.5625) 


(0.5300, 0.5650) 


First BL 


^max 




2.12504E-4 


2.31129E-4 


2.3261E-4 




Location (x, 


y) 


(0.0781, 0.0781) 


(0.0859, 0.0781) 


(0.0833, 0.0783) 


First BR 


^max 




1.67313E-3 


1.75102E-3 


1.7281E-3 




Location (x, 


y) 


(0.8672, 0.1094) 


(0.8594,0.1094) 


(0.8633, 0.1117) 


Second BR 






-4.815059E-8 


-9.31929E-8 


5.4962E-8 




Location (x, 


y) 


(0.9922, 0.0078) 


(0.9922, 0.0078) 


(0.9917, 0.0067) 



Table 21: 


Streamfunction values 


*mm, ^max and locations 


of the primary and secondary 


vortices 








Number 


property 


Re=3200 


Re=3200 






Multistep scheme with index 


m = 2 Ghia et al. [27] 


Primary 




-0.109962 


-0.120377 




Location, x, y 


(0.5156,0.5391) 


(0.5165,0.5469) 


First T 


^max 


5.759079E-4 


7.27682E-4 




Location (x, y) 


(0.0469,0.8984) 


(0.0547,0.8984) 


First BL 


^max 


1.09512E-3 


9.7823E-4 




Location (x, y) 


(0.0781,0.1250) 


(0.0859,0.1094) 


First BR 


^max 


2.70425E-3 


3.13955E-3 




Location (x, y) 


(0.8281,0.0859) 


(0.8125,0.0859) 


Second BL 




-1.04040E-8 


-6.33001E-8 




Location (x, y) 


(0.0078,0.0078) 


(0.0078,0.0078) 


Second BR 




-1.36461E-7 


-2.51648E-7 




Location (x, y) 


(0.9844,0.0078) 


(0.9844,0.0078) 
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-0.4 -0.2 0.2 0.4 0.6 0.8 1 -0.5 0.5 1 



(a) Re = 1000 (b) Re = 3200 

Figure 1: Velocity (ux) profiles along the vertical line passing through the geometric center of 
the cavity. Black solid lines: (a) single-step scheme, (b) multistep scheme with index m = 2; 
Blue circle lines: Ghia et al. |27| 




0.2 0.4 0.6 0.8 1 ■ 0.2 0.4 0.6 0.8 1 

X y 



{a)Re = 1000 {h)Re = 3200 

Figure 2: Velocity (uy) profiles along the horizontal line passing through the geometric center 
of the cavity. Black solid lines: (a) single-step scheme, (b) multistep scheme with index m = 2; 
Blue circle lines: Ghia et al. |27j 



time step size, which can be essentially improved by an explicit multistep scheme with smaller 
time step sizes so that the resulting method behaves like an unconditionally stable method. The 
diffusion problem involved at each iteration is always self-adjoint and coercive so that it can 
be solved efficiently using many existing optimal preconditioned iterative solvers. The optimal 
convergence orders have been confirmed by several numerical examples with smooth solutions. 
The schemes are then extended for the Navier-Stokes equations, where the nonlinearity is re- 
solved by a linear explicit multistep scheme at the convection step, while only a generalized 
Stokes problem is needed to solve at the diffusion step and the major stiffness matrix stays 
invariant in the time marching process. Numerical simulations are presented to demonstrate 
the stability, convergence and performance of the single-step and multistep variants of the new 
schemes. The effectiveness and robustness of the new schemes are finally well demonstrated by 
the benchmark lid-driven cavity fiow problem. The newly proposed schemes are all free from 
tuning some stabilization parameters as the most existing schemes require for the convection- 
dominated diffusion problems. Finally we note that the proposed fully discrete schemes are only 
first order in time, and we are currently investigating the potential schemes which have second 
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0.2 0.4 0.6 0.8 1 -5 5 10 15 20 

K vorticity 



Figure 3: Vorticity values along the vertical line x = 0.5 (left) and the horizontal line y = 0.5 
(right) passing through the geometric center of the cavity with Re = 1000. Black solid lines: 
single-step scheme; Blue circle lines: Botella et al. [29] 

order temporal accuracy. 
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